
## functions
plot.coef <- function(modResults, data, vars) {

modelcoef=lapply(modResults, 
                 function(x) FUN=summary(x)$coefficients[2:length(x$coefficients), 1])
modelse = lapply(modResults, 
                 function(x) FUN=summary(x)$coefficients[2:length(x$coefficients), 2])

varnames =lapply(modResults, function(x) FUN = names(x$coefficients)[-1])

dfplot <- list()

for (i in 1:length(modelcoef)){
       ad = 0.5 / modelse[[i]]
       modelcoef[[i]] = modelcoef[[i]]* ad
       modelse[[i]] = modelse[[i]]*ad 
       
       ylo95 =modelcoef[[i]] - qt(.975, nrow(data))*(modelse[[i]])
       yhi95 =modelcoef[[i]] + qt(.975, nrow(data))*(modelse[[i]])
       ylo90 =modelcoef[[i]] - qt(.95, nrow(data))*(modelse[[i]])
       yhi90 =modelcoef[[i]] + qt(.95, nrow(data))*(modelse[[i]])
       dfplot[[i]] = data.frame(varnames[[i]], modelcoef[[i]], modelse[[i]], ylo95, yhi95,
                                ylo90, yhi90,
                                Model = paste('Model',seq_along(modelcoef)[i]))
}
coefficientdata = do.call(rbind,dfplot)
colnames(coefficientdata) <- c("Variables", "Coefficients", "S.E",
                               "Low95CI", "High95CI", 
                               "Low90CI", "High90CI",
                               "Model.Name")

df <- coefficientdata %>%
                     dplyr::mutate(Variables = as.character(Variables)) %>%
                     dplyr::arrange(Variables)

df <-  df[grep(paste(vars, collapse = "|"), df$Variables),]
       

p = ggplot(df, aes(colour =Model.Name)) + 
       geom_hline(yintercept = 0, lty = 2) +
       geom_linerange(aes(x = Variables, ymin = Low90CI,
                          ymax = High90CI),
                      lwd = 1, position = position_dodge(width = 1/2)) +
       geom_pointrange(aes(x = Variables, y = Coefficients, ymin = Low95CI,
                           ymax = High95CI,shape = Model.Name),
                       lwd = 0.5, position = position_dodge(width = 1/2),
                       fill = "WHITE") +
       coord_flip() + 
       theme_bw() + xlab("") + ylab("Rescaled Coefficient") + 
       theme(legend.position="bottom",
             legend.title=element_blank(),
             axis.text = element_text(size=14),
             text = element_text(size=14),
             plot.title = element_text(hjust = .5, size = 16, face = "bold"))
return(p)   
}  


###Coef plot for table 2
plot_interEffect = function(ModelResults, n.sim = 1000, data, varname1, varname2,
                            val1, val2, intervals, label = c("lab1", "lab2"), xlabs, ylabs, facet = F){
       #get a sim objective
       require(arm)
      # library(multiwayvcov)
       #library(lmtest)
       library(ggplot2)

       set.seed(12345)
       sim <- mvrnorm(n= n.sim, coef(ModelResults), vcov(ModelResults))
       
       ##set simulation
       varname2_val = seq(val1, val2, by =intervals)
       df_no <- array(NA, c(length(varname2_val), 5))
       df_yes <- array(NA, c(length(varname2_val), 5))
       
       for (i in 1:length(varname2_val)){
              X1 <- model.matrix(ModelResults)
              X2 <- model.matrix(ModelResults)
              X1[, varname1] = 0
              X1[, varname2] = varname2_val[i]
              X1[, paste(varname1, varname2, sep = ":")] = 0*varname2_val[i]
              
              X2[,varname1] = 1
              X2[,varname2] = varname2_val[i]
              X2[, paste(varname1, varname2, sep = ":")] = 1*varname2_val[i]
              
              fd_no = apply(apply(X1, 1, function (x) plogis(sim %*% x)), 1, mean)
              
              fd_yes = apply(apply(X2, 1, function (x) plogis(sim %*% x)), 1, mean)
              
              df_no[i, 1] <- varname2_val[i]
              df_no[i, 2] <- mean(fd_no)
              df_no[i, 3:4] <- quantile(fd_no, probs = c(.05,.95))
              df_no[i, 5] <- 0
              df_yes[i, 1] <- varname2_val[i]
              df_yes[i, 2] <- mean(fd_yes)
              df_yes[i, 3:4] <- quantile(fd_yes, probs = c(.05,.95))
              df_yes[i, 5] <- 1
       }
       
       df_plot <- rbind(df_no, df_yes)
       colnames(df_plot) <- c("X", "mean", "lo", "hi", "type")
       df_p <- as.data.frame(df_plot)
       df_p$type <- as.factor(df_p$type)
       levels(df_p$type) <- label
       
       if(facet == F){
              p =   ggplot(df_p, aes(x=X, y=mean, group = type, color = type, linetype=type)) +
                     theme_gray() +
                     geom_ribbon(aes(ymin = lo, ymax = hi), alpha=.2, color = NA) +
                     geom_line(aes(y = mean, x = X), size = 0.5, alpha = 0.75) +
                     ylab(ylabs) + theme_bw() +
                     scale_color_manual( "", labels= label,
                                         values=c("#2c7bb6","#d7191c")
                     ) +
                     scale_linetype_manual("",
                                           labels=label,
                                           values=c(1,2)) +
                     scale_x_continuous(name = xlabs, breaks = seq(val1, val2, by = 2)) +
                     theme(axis.title.y = element_text(vjust=1,size=12),
                          # axis.title.x = element_blank(),
                           legend.position="bottom",
                           legend.margin=margin(0,0,0,0),
                           legend.text=element_text(size=12),
                           axis.text.y=element_text(size=12),
                           axis.title=element_text(size=14),
                           axis.text.x=element_text(size=12),
                           plot.margin=unit(c(0.1,0.1,0.1,0.1),"cm"))
              return(p)}
       else { p =   ggplot(df_p, aes(x=X, y=mean)) + theme_gray() +
              geom_ribbon(aes(x = X, ymin = lo, ymax = hi), alpha=.4) +
              stat_smooth(color = "#2c7bb6") + facet_wrap(~type, scales = "free_y") +
              scale_x_continuous(name = xlabs, breaks = seq(val1, val2, by = 2)) +
              ylab(ylabs) + theme_bw() +
              theme(axis.title.y = element_text(margin = margin(1,1,1,1)),
                    axis.text = element_text(size=14),
                    axis.title=element_text(size=14),
                    strip.text = element_text(size=15))
       return(p)
       }
       
}

###define a theme in ggplot2
theme_nothing <- function(base_size = 12, legend = FALSE){
   if (legend) {
      return(theme(line = element_blank(),
                   rect = element_blank(), #defien the margin line
                   axis.text = element_blank(), 
                   axis.title = element_blank(), 
                   panel.background = element_blank(),
                   panel.grid.major = element_blank(), 
                   panel.grid.minor = element_blank(),
                   axis.ticks.length = unit(0,"cm"),
                   panel.spacing = unit(0, "lines"),
                   plot.margin = unit(c(0,  0, 0, 0), "lines")))
   }
   else {
      return(theme(line = element_blank(), 
                   rect = element_blank(), 
                   text = element_blank(),
                   axis.ticks.length = unit(0,"cm"), 
                   legend.position = "none",
                   panel.spacing = unit(0, "lines"),
                   plot.margin = unit(c(0, 0, 0, 0), "lines")))
   }
}